Packing geometry and statistics of force networks in granular media 
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The relation between packing geometry and force network statistics is studied for granular media. 
Based on simulations of two-dimensional packings of Hertzian spheres, we develop a geometrical 
framework relating the distribution of interparticle forces P{f) to the weight distribution V{w), 
which is measured in experiments. We apply this framework to reinterpret recent experimental data 
on strongly deformed packings, and suggest that the observed changes of 'P{w) are dominated by 
changes in contact network while P{f) remains relatively unaltered. We furthermore investigate 
the role of packing disorder in the context of the g-model, and address the question of how force 
fluctuations build up as a function of the distance beneath the top surface. 



PACS numbers: 45.70.-n, 45.70.Cc, 46.65.+g, 05.40.-a 



I. INTRODUCTION 

Inside a granular material forces are distributed very 
inhomogeneously: a small number of particles carries a 
large fraction of the internal forces [H . These large fluctu- 
ations are reflected in the force probability density func- 
tions, which typically decay exponentially 0, S S 0|- 
The behavior for small forces is not as well understood as 
the generic exponential tail: the g-model appears to pre- 
dict a vanishing probability density for small forces [3, 
whereas experiments and simulations clearly show that 
this probability remains non-zero |^ |^ Q . The charac- 
terization and understanding of this probability remains 
a challenge, especially since the force distribution is be- 
lieved to play an important role for the dynamical arrest 
or ''jamming" of granular and other disordered materi- 
als In particular, the force distribution has been ob- 
served to develop a small peak (around the average value) 
in simulations of supercooled liquids, foams and granular 
matter undergoing a jamming transition However, 
there is still no microscopic understanding how this effect 
relates to the properties of the force network. 

This paper is a full exposition and expansion of a new 
approach, which was briefly outlined in We will un- 
ravel the effect of the local contact geometry on the dis- 
tributions of interparticle force F and effective particle 
weight W; the weight is defined as the sum of the ver- 
tical components of all downward pointing forces on a 
particle - see Fig. 1. While the distribution of forces F is 
the primary object one ultimately wishes to characterize, 
it is difficult to access experimentally. Experiments with 
photoelastic materials are able to depict the spatial struc- 
ture of bulk forces in 2D, but their precision to resolve 
individual contact forces is limited Only recently, 
there have been first reports of 3D bulk measurements 
on forces in compressed emulsions p^. Most quantita- 
tive information on the force probability distribution is 
at present only accessible through measurements of the 
particle-wall forces from imprints on carbon paper 2] or 
by force sensors Q. Each particle- wall force has to bal- 
ance all interparticle forces that are exerted on the corre- 
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FIG. 1: (a) Detail of a typical packing in our simulations; 
the height h denotes the distance from the bottom. The force 
network is represented by the black lines whose thickness is 
proportional to the force-magnitude, (b) Definition of inter- 
particle forces F and weight W , for a frictionless particle with 
nc — 2. 



sponding particle from above, see Fig.^ This means that 
experiments essentially measure a combination of forces 
that we refer to as the weights of the bottom particles. 
For simplicitly, we will focus on frictionless spheres for 
which these weights are defined as 



mjg 



(1) 



Here mj denotes mass, g denotes gravity, Fij are the 
interparticle forces and ric is the number of particle ex- 
erting a force on particle j from above; the sum runs over 
all these forces. So, to relate the experimental results to 
the bulk force distributions, one has to understand the 
relation between weights and forces. 

In this paper we will show how the local packing geome- 
try plays the crucial role in the relation between the force 
distributions P{f) and the weight distributions V{w) (we 
define f = F/ (F) and w — W/ {W) as the appropriately 
rescaled forces and weights). Our central point is that 
while the distribution of / is robust, the distribution of 
w is profoundly influenced by the contact geometry, in 
particular by the number of downward pointing contact 
forces ric . In simulations of Hertzian sphere packings we 
will find that V boundary {w) is different from Vbuikiw), due 
to the rather special packing geometry near a boundary. 
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However, for many (but not all) experimentally relevant 
situations, the special packing geometry near a boundary 
makes Vboundaryiw) rather close, but not equal, to the 
bulk P{f). This fortunate but non-trivial coincidence 
can be understood easily within our framework. We will, 
however, also provide two examples where Vboundaryiw) 
and bulk P{f) are significantly different. 

Additional motivation for studying the relation be- 
tween forces, weights and geometry comes from the q- 
model Once the distinction between forces and 

weights has been made, one notices that the g-model is a 
lattice model in which weights are randomly redistributed 
over a fixed number of supporting grains. The g-model 
displays a weight distribution that is qualitatively dif- 
ferent from both experimentally observed weight distri- 
butions, or numerically obtained force distributions. We 
will show that this is due to the fixed connectedness of 
the q-model. Realistic V{w) can be obtained if we allow 
for the connectivity to vary within the g-model, e.g. by 
introducing random connectivity. 

Our work then serves three purposes. First of all, 
it helps to interpret data obtained by measurements of 
particle-wall forces: this paper includes a section where 
we explicitly apply our framework to recent experimen- 
tal data of highly compressed packings ■ Secondly, it 
shows how the simple g-model can be extended to ob- 
tain very realistic weight distributions for both regular 
and irregular packings. Since the model is known to give 
incorrect predictions on spatial propagation [T^ . our in- 
tention is not to fine-tune the model and its parameters, 
but rather to indicate how the contact geometry is essen- 
tial to describe force and weight fluctuations in more re- 
alistic packings. Thirdly, we address the question of how 
force fluctuations build up as a function of the distance 
beneath the top surface, providing another fundamental 
test for theoretical models. 

The paper is organized as follows. In Sec. we first 
explain our numerical method and then discuss the force 
distributions observed in amorphous packings: it turns 
out that P{f) is rather insensitive to the packing geom- 
etry. We then show in Sec. IIIII that the weight distri- 
butions V(w), on the other hand, are very sensitive to 
the packing geometry. Using simple phase space con- 
siderations, we relate V{w) to P{f) for a given geome- 
try. This provides a recipe how to reconstruct the bulk 
P{f) from the experimental data, and in Sec. IIVI wc ex- 
plicitly apply this to recent experimental data on highly 
compressed packings particular, our analysis 

strongly suggests that P{f) is essentially unaffected by 
the tremendous deformations encountered in the experi- 
ments. We then indicate some limitations of our frame- 
work in Sec.0 where we address subtle packing problems 
like the effect of gravity. In Sec. I VII we investigate to what 
extent the q-model can describe the results of the numer- 
ical packings of Hertzian spheres: we derive a surprising 
exact result for the bond quantities qw, and we investi- 
gate the role of disorder in the packing geometry. Finally, 
we address the top-down relaxation of force fluctuations 



in Sec. IVIII We find no evidence in the Hertzian sphere 
packings for the power-law relaxation predicted by the 
g-model, indicating that the model is not able to capture 
this spatial aspect of the force network. The paper ends 
with a discussion. 



II. STATISTICS OF INTERPARTICLE FORCES 

In this section we study the distribution of interpar- 
ticle forces via simulations of 2D packings of friction- 
less spheres. After introducing our numerical method in 
Sec. Ill AI we discuss the similarities between P{f ) in the 
bulk and near the boundary fSec. Ill B)| . We also study 
the angular distribution and the probability distribution 
of the z components of the contact forces in Sec. Ill CI 
and close with a brief summary of results in section FlIDI 



A. Numerical method and parameters 

Our two-dimensional packings consist of frictionless 
spheres (3D) under gravity. The packings are created 
from molecular dynamics simulations of spheres that in- 
teract through normal Hertzian forces, where F oc 6^^'^ 
and S denotes the overlap distance 0|. Since Hertz's law 
for 2D disks is linear in 5, we use 3D spheres. These par- 
ticles reside in a container that is 24 particle diameters 
wide, with periodic boundary conditions in the horizon- 
tal direction. The bottom support is rigid and also has 
a frictionless Hertzian interaction with the particles. We 
construct our stationary packings by letting the particles 
relax from a gas-like state by introducing a dissipative 
force that acts whenever the overlap distance is non-zero. 
In this paper we use two different polydispersities: the 
radii r are drawn from a flat distribution between ei- 
ther 0.49 < r < 0.51 or 0.4 < r < 0.6. The masses 
are proportional to the radii cubed. In the former case 
of almost monodisperse particles, the particles tend to 
crystallize into a triangular lattice fSec. IIVX)l . whereas 
the more polydisperse particles lead to amorphous pack- 
ings such as shown in Fig. QJi. This allows us to study 
how the packing geometry affects the force network. The 
results shown in this paper are obtained with particles 
that deform 0.1% under their own weight. Simulations 
of harder particles (deformation 0.01%) gave similar re- 
sults as those shown here . 

The various data were obtained from 1100 realizations 
containing 1180 particles each. We study the force and 
weight distributions at various heights h. To do so, we 
divide each packing into horizontal slices of one particle 
diameter thickness, and rescale all forces and weights in 
each layer to the corresponding average (absolute) val- 
ues. The rescaled interparticle forces and weights will be 
denoted by / and w respectively, with distributions P{f) 
and V^w). 
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FIG. 2: (a) P{f) for amorphous packing in the bulk (open cir- 
cles) and for the layer-to-layer forces near the bottom (dots); 
the inset shows P{f) for bulk forces on a log-lin scale. Note 
that the force distributions are very similar, except for a small 
difference for small /. (b) Detail of a typical packing near 
the bottom showing layer-to-layer forces (black lines) and the 
intralayer forces (white lines) near the bottom. It is clear 
that the layer-to-layer forces are dominant in determining the 
weights w of the bottom particles. The numbers show the val- 
ues of He, the number of (layer-to-layer) forces that contribute 
to these weights. 



B. Absolute values of /: P(f) 

We first analyze the statistics of the absolute values 
/ = I/I J whose probability density function P(/) is 
usually referred to as the distribution of (interparticle) 
forces; our main finding will be that P{f) in bulk and 
near the boundary are very similar. In Fig. we show 
P{f) as measured in the bulk of the amorphous pack- 
ings (particle radii between 0.4 < r < 0.6). At different 
heights between 10 < ft. < 30, P{f) was not observed to 
change; the open circles represent an average over these 
various heights. Even very close to the bottom support, 
we find that P{f) remains almost unchanged: the dot- 
ted dataset has been obtained from the forces between 
the bottom particles and the particles in the layer above. 
We refer to these forces as layer-to-layer forces near the 
bottom - see Fig. I^Jj). So, although the bottom wall lo- 
cally alters the packing geometry, the shape of P{f) is 
essentially unaffected. 

As can be seen from the inset of Fig. |21 the probability 
density decays slightly faster than exponentially. This 
is consistent with simulations by Makse et al. [I3 who 
found that P{f) crosses over to a Gaussian for large par- 
ticle deformations; we have used rather 'soft' particles 
in our simulations for which deformations are relatively 
large, i.e. up to 2%. We come back to the effect of de- 
formation in experiments in Sec. lIVBl For small forces, 
P{f) approaches a finite value. The small peak around 
/ = 0.7 for bulk forces becomes a plateau for the layer- 
to-layer forces near the bottom; it is intruiging to note 
that this change is reminiscent of what is proposed as an 
identification of the jamming transition 
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FIG. 3; Scatter plot of {fij,ipij) for (a) the bulk forces, and 
(b) the layer-to-layer forces near the bottom in the amorphous 
packings. 



C. Orientations of / and P'{fz) 

After studying the absolute values of /y , let us in- 
vestigate the orientations of the interparticle forces. We 
therefore define Lpij as the angle between fij and the hori- 
zontal axis. In Fig.|3^ we show the scatterplot of {fij , fij ) 
in the bulk: the angles are uniformly distributed and in- 
dependent of the absolute value of /. So, the packings 
are highly disordered away from the bottom. Near the 
boundary, however, this isotropy is broken strongly. The 
presence of the bottom wall aligns the bottom particles 
and as a consequence their interparticle forces become al- 
most purely horizontal, see Fig. \Bp. It is clear that near 
the bottom the interparticle forces naturally divide up 
into these almost horizontal intralayer forces, and layer- 
to-layer forces connecting bottom particles with those in 
the layer above. The orientations of these layer-to-layer 
forces are indeed concentrated around 7r/3 and 27r/3, as 
can be seen from Figs.l^t and|3|D. 

Since the particle weights are derived from the z- 

components of the forces, fz = i^fij^ ■> now investi- 
gate their distribution P'{fz). The bottom-induced ori- 
entational order discussed above is reflected in the statis- 
tics of the fz ■ According to Fig. ^ there is a substan- 
tial difference between P'{fz) in the bulk (open circles) 
and P'ifz) for the layer-to-layer forces near the bottom 
(dots). This difference can be understood as follows. As- 
suming that the tpij are indeed uncorrelated to the fij, 
we can write 



dip^f) dfP[f)5{fz^fsHv)). 

JQ 



(2) 

where ^{(p) is the angle distribution, and P{f) is the 
distribution of the absolute values |/| of Fig. |21 Note 
that {fz) < 1. For the layer-to-layer forces near the bot- 
tom, we have seen from the scatter plot that the values 
of sin{if) are concentrated around ~ 0.866. In the 
approximation that the distribution of sin{ip) is sharply 
peaked, the shape of P'{fz) equals that of P{f) (up to a 
scale factor). This is indeed conflrmed by direct compar- 
ison of the dotted datasets of Figs. |21 and ^ 

In the bulk, we have seen that the packing geome- 
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FIG. 4: P'ifz) in the bulk (open circles) and for the layer-to- 
layer forces (dots). The solid line was obtained by numerical 
integration of Eq. 0. Inset shows P'ifz) versus log fz, con- 
firming the logarithmic divergence for small /z. 



granular packing, and not the interparticle (bulk) forces 
that were discussed in the previous section. Since these 
particle- wall forces are essentially equal to the weights of 
the bottom particles, it is important to understand the 
relation between the weight distribution 'P{w) and the 
distribution of interparticle forces P{f )- In the first part 
of this section we develop a simple geometrical frame- 
work to understand this relation, based on phase space 
considerations. We then show that this explains, to a 
large extent, the weight distributions V{w) as measured 
in our simulations of Hertzian spheres. In particular, we 
observe substantial differences between weight distribu- 
tions for different packing geometries. 



try is isotropic. A consequence of this isotropy is that 
the probability density function of the horizontal com- 
ponents P'ifx) is identical to P'ifz) (not shown here). 
Again, one can use Eq. (|2Jl to understand the shape of 
P'ifz)- Taking a uniform angle distribution <i>(<p) = I/tt, 
we obtain (Appendix 

^'(/^)-^ r'^/^==f (3) 

^ Jf. Vf - fz 

Numerical integration of this equation with P(/) from 
Fig. 12 yields the solid line in Fig. 01 which closely cor- 
responds to the P'ifz) as measured in the bulk (open 
circles) . In Appendix we show that the integral of 
Eq. (|3J) is weakly divergent for small fz- 

P'ifz)^--PiQ)Hfz) + Oil) . (4) 

TT 

The inset of Fig. 0] shows that our data for P'ifz) is 
indeed consistent with this logarithmic divergence. 

D. Pif)- summary 

Let us briefly summarize the results of this section. 
The geometrical constraint imposed by the bottom wall 
locally induces a packing geometry which is different from 
the bulk geometry. Whereas this is strongly reflected in 
the orientations of the fij, the distribution of the ab- 
solute values Pif) is very robust. The probabilities for 
the components of the fij can be obtained with great 
precision, including the logarithmic divergence, by the 
transformation of Eq. 

III. PACKING GEOMETRY AND WEIGHT 
DISTRIBUTIONS P(to) 

In this section, we demonstrate that the local packing 
geometry has a dramatic effect on the weight distribution 
of P(w). As stated in the introduction, experiments can 
only measure the particle- wall forces at the boundary of a 



A. Geometrical framework: decomposition of Viw) 
according to number of contacts ric from above 

If we interpret Eq. as a transformation of stochas- 
tic variables, it is possible to relate the corresponding 
probability density functions as 

VnAW) - / diFi)z--- diFnJz 
JO JO 

X p[iFi)zr--.iFnJz) s(w^p^iF,)}j . 

(5) 

Here, we have neglected the term mg, since mg/ (W) <C 1 
far below the top surface of the packing. The number of 
forces over which we integrate differs from grain to grain, 
and it turns out to be crucial to label the weight distri- 
bution in Eq. ((SJ, Vn^iW), according to this number Uc- 
This can be seen as follows. The (5-function constrains the 
integral on a (n^ — 1) dimensional hyperplane of the total 
phase space, and the "area" of this hyperplane scales as 
T4^"<:-i. We thus anticipate the following scaling behav- 
ior for small weights: 

Vn^W) (xW'^"-'^ for w^O, (6) 

provided that the joint probability density approaches 
a finite value when all iFi)z 0. Such scaling is also 
implicit in the g-model Q , although there ric > 2 so that 
7^(0) = 0. The particles that do not feel a force from 
above, ric = 0, give a (5-like contribution at W = mg; for 
deep layers this occurs for mg/ {W) ^1. In a disordered 
packing, the number of particles that exert a force from 
above can vary from grain to grain. The total weight 
distribution ViW), therefore, is a superposition of the 

'PnAW). 

ViW)^Y.Pr^^^r.M) , (7) 

Wo 

where pn^ is the fraction of particles with Uc contacts 
from above. This means that the small weight behavior 
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of ■p(VF) depends very much on the fractions and 
thus on the local packing geometry, via Eqs. © and iQ. 

The steepness of the tail of the total weight distri- 
bution depends strongly on the as well. To ex- 
plain this, let us assume that all vertical forces Fz con- 
tributing to the weight are uncorrelated. We consider 
P'Uz) oc e-"'!', i.e. P'(F^) oc e'"'^-/'^^'^ for large 
forces. It follows from Eq. (O that the weight distri- 
bution takes over this same exponent a/{Fz), so that 
VnAW) OC e-"^/<^->. However, the Pn,(W^)'s are not 
properly normalized: {W)n^ — {Fz) n^, since each of the 
Fz gives an average contribution (Fz). This yields a to- 
tal average weight (W) = {Fz)J2„, Pn^^c = {Fz){nc). 
In order to compare with experimental and theoretical 
results we have to rescale the weights so that (w) = 1, 
yielding the following large weight behavior: 

Viw) oc e"^*" with 7 = a (n^) • (8) 

This simple calculation shows that, for a given value of 
a, the steepness of the tail of the experimentally mea- 
sured weight distribution is very sensitive to the local 
packing geometry. This is a direct consequence of keep- 
ing (w) fixed to unity: a decrease of probability for small 
weights must lead to a steeper tail for large weights in 
order to leave the average weight unaltered. Note that 
this general argument is not restricted to uncorrelated 
Fz or exponential tails. A generalization to other than 
exponential tails is given in Appendix IbI 

So, we have advanced a simple picture, in which the 
shape of Viw) depends strongly on the local packing ge- 
ometry via the fractions The small force behavior 
follows from Eqs. © and Q), whereas Eq. ||SJ) relates to 
a good approximation the exponential tails of P' [fz) and 
V{w). The object one ultimately wishes to characterize 
is of course the force distribution P{f). Since close to 
the boundary P{f) and P'{fz) are identical up to a scal- 
ing factor {fz) (Sec. Ill the above equations allow to 
trace the features of the force distribution from experi- 
mental measurements. Along this line, we analyze recent 
experimental data in Sec. IIV Bl 

B. Viw) in Hertzian sphere packings 

We now discuss the weight distributions observed in 
the Hertzian sphere packings, and interpret the results 
within the framework developed above. Figure shows 
that in the amorphous packing V{w) in the bulk (open 
circles) is significantly different from 7^ (it;) of the bot- 
tom particles (dots). The probability for small weights 
is much larger at the bottom, and the decay for large 
weights is not as steep as for the bulk particles. Fur- 
thermore, the transition from bottom to bulk behavior is 
remarkably sharp: in the slice 2 < ft, < 3 (full curve), the 
weight distribution is already bulk-like. 

Using the concepts developed in the preceding para- 
graphs, we now show how this change in 'P{w) can 
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FIG. 5: (a) V{w) in the bulk (open circles) and at the bottom 
(dots) in amorphous packings. At 2 < /i < 3, V{vj) is already 
bulk-like (solid line). Decomposition of V(w) accord- 

ing to Eq. Q (V) in the bulk (open circles) and (c) a,t the 
bottom (dots). The measured bulk values for the fractions 
{po,pi,p2,P3} in Eq. ^ are {0.01,0.11,0.52,0.36}, and the 
bottom values are {0.08, 0.46, 0.44, 0.02}; as explained in [3, 
we excluded the intralayer (almost horizontal) forces at the 
bottom when determining n^. 



be explained by a change in the local packing geom- 
etry. Consider the typical bottom configuration of 
Fig. The intralayer forces (white lines) are almost 
purely horizontal and hence do not contribute to the 
weights. This reduces the effective values of Uc, lead- 
ing to the following fractions for the bottom particles: 
{po,Pi,P2,P3} = {0.08,0.46,0.44,0.02}, where we did 
not count the intralayer forces for determining the val- 
ues of He 01 . In the bulk, these fractions are different, 
namely {po, Pi, P2, Ps} = {0.01,0.11,0.52,0.36}. Accord- 
ing to Eq. (UJ, these differences between the in the 
bulk and at the bottom should lead to a substantially 
different V{w). Figs.[SlD,c exphcitly shows the decompo- 
sition into the VnA^)- Indeed, one observes the scaling 
behavior for small w proposed in Eq. ((HJ. Moreover, the 
various Vn^{w) are essentially the same at the bottom 
and in the bulk: a direct comparison is given in Fig. El 
where we rescaled the average values to unity. There 
is only a small difference in the Vi{w) due to the fact 
that bottom particles with ric = 1 are typically smaller 
than average (Fig.^t). For these particles, the intralayer 
forces will add a small contribution to the weights, en- 
hancing Vi{w) for small w at the expense of 'Pi(O). The 
same argument holds for Vo{w), whose (5-like shape ap- 
pears a bit broadened in Fig. However, it is clear 
that the differences between Viw) in the bulk and at the 
bottom are mainly due to a change in contact geometry. 

Finally, let us remark that the good agreement between 
Pbuikif) and Vboundaryiw) for w > 0.3 is fortuitous and 
due to the relatively large fraction of bottom particles 
with Tie = 1. We will argue below that this is also the 
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FIG. 6: Direct comparison of (a) Viiyo) and (b) V2(w) for 
bulk (open circles) and bottom particles (dots). All distribu- 
tions are scaled such that (w) = 1. 



case in many (but not all) carbon paper experiments. 



C. Summarizing the simple picture 

Our simple framework as developed in the sections 
above can be summarized as follows: The geometry of 
the contact network has a strong effect on P{w), while 
P(/) is very robust. The weight distribution for par- 
ticles with a given Uc, Vn^{w), is robust and behaves 
as w"''^^ for small w. V{w) can be decomposed as 
T^iw) = J2n Pric'Pn^i'w), whcrc are the fractions of 
particles that have ric = 0, 1, 2, . . . "up" contacts. Differ- 
ences of between boundary particles and bulk parti- 
cles explain the different 7-'(w)'s for these cases. When 
po and pi are large, the total weight distributions 7^(it;) 
exhibits a plateau at small weights and a slow decay at 
large weights; when p2 and p^ become large, V{w) be- 
comes sharply peaked. In this way, 7'(m;)'s small weight 
behavior as well as its exponential decay rate for large 
weights reflect the packing geometry. 



IV. MANIPULATING THE GEOMETRY: 
EXPERIMENTAL RELEVANCE 

So far we have focused on the role of the bottom bound- 
ary for disordered packings of frictionless particles. In 
this section we provide explicit examples of other types 
of packing geometries and their effect on 'P{w). We first 
discuss our simulations of weakly polydisperse particles, 
which give rise to rather crystalline packings - see Fig.Et^. 
We then apply the geometrical framework derived in the 
previous section to experimental (carbon paper) data by 
Erikson et al. TTj of highly deformed packings of soft 
rubber particles. Their results have a natural interpre- 
tation within our framework and form a nice illustra- 
tion of how the number of contact affects the weight dis- 
tribution. Both the simulations of crystalline packings 
and the experiments on deformed packings are examples 
where the experimentally accessible Vboundaryiw) is sig- 
nificantly different from P{f) in the bulk; we discuss why 
in many other carbon paper experiments Vhoundary{w) is 
probably very similar to the real P{f). 
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FIG. 7: (a) Weakly polydisperse particles (radii between 
0.49 < r < 0.51) spontaneously crystallize into a hexago- 
nal packing, (b) The corresponding P{f) is indistinguishable 
from the force distributions in amorphous packings, (c) The 
weight distributions V{w) in the bulk (open circles) and at 
the bottom (dots) are dominated by particles with — 2. 



A. Crystalline versus disordered frictionless 
packings 

We now present the results of the more or less crys- 
talline packings, obtained from simulations with parti- 
cle radii between 0.49 < r < 0.51. Firstly, the force 
distribution P(/) shown in Fig. \7}p is indistinguishable 
from the force distributions in the amorphous packings 
(compare with Fig.|2t). So despite the order in particle 
positions, there are still large fluctuations in the force 
network. There is of course some disorder in the "con- 
tact network" since not all particles are in contact with 
their six neighbors (Fig. [7^). It is nevertheless surprising 
that for this very different contact geometry, the force 
fluctuations are characterized by the same probability 
distribution as was observed for highly disordered pack- 
ings. This strongly suggests that P{f) is a very robust 
quantity and independent of the packing geometry. 

The weight distribution V(w), on the other hand, is 
very sensitive for the geometry. In a perfect triangular 
packing all particles would have ric = 2; in our simula- 
tions we find that p2 = 0.9 and pi =0.1 due to lattice im- 
perfections. From our geometrical framework we expect 
that the shape of the weight distribution is dominated by 
7^2 (w). Fig. [7|: shows that this is indeed the case - e.g. 
compare with Fig. |^. 

In an earlier paper Q , we reported how one can break 
the regular packing geometry by using curved bound- 
aries. This led to a dramatic change in V{w) that again 
could be understood from a change in the . 



B. Experiments on strongly deformed particles 

We now demonstrate how the strategy to decompose 
the weight distributions according to ric can be applied to 



7 



experiments measuring V{w) at the boundary of a gran- 
ular material. This is best illustrated by recent carbon 
paper experiments by the Chicago group on soft rubber 
beads, in particular Fig. 3 of Ref. [ill, which the effect 
of particle deformations was investigated. The raw data 
of these experiments were kindly made available by the 
authors, allowing us to perform the analysis presented 
below. 

The experimental results of Fig. 3 of Ref. display 
three trends as the compression is increased: 

• The (5-like peak at w = decreases, 

• lmijjuiQV{w) decreases, 

• The exponential tail becomes steeper. 

These behaviors emerge naturally when considering the 
role of the fractions /5„^. The first trend arises from a 
decrease in po, since only particles with Uc — give a 6- 
like contribution to Viw). The second trend comes from 
a decrease in pi. from Eqs. © and Q it is clear that 
\\m.u,ioV{w) — piVi{w). The changes in V{w) can thus 
be understood from an increasing number of contacts, 
which is what one would expect for a compressed system 
|15| . The fractions p2 and p^ will increase at the expense 
of pq and pi. Also the third trend, the steepening of 
the exponential tail, is directly related to the increase in 
(ric) via Eq. ||SJ). However, Eqs. I©-© allow to further 
quantify this change in contact geometry from the exper- 
imental data. The value of piViiO) can be read off from 
the plots, after subtracting the 5-like data points, since 
Pi^i(O) = \im^^QV{w). The value of pq is obtained by 
the height of the S-peak times the bin-width. Using the 
raw experimental data, we obtained the figures given in 
the first colomn of Tabled where we took 7^i(0) = 0.5 
|17| . Unfortunately, the values of p2 and ps can not be 
determined directly from the data. 

An intriguing issue is that numerical simulations by 
Makse et al. |13| indicate that P{f) crosses over to a 
Gaussian for large particle deformations. This contra- 
dicts the experimental data for which one observes an 
exponential tail even though particle deformations are up 
to no less than 45% 0. Moreover, we speculate below 
that the steepening of the tails is only due to changes in 
the and that the bulk force distributions P{f) actu- 
ally remain unaffected by the particle deformations. The 
way to test this scenario is to examine whether the ex- 
ponential decay constant of P{f) oc e~°'^ remains fixed, 
even though the steepness oiV{w) oc e"''™ increases. We 
use Eq. ^ to determine the value of a = ^/{ric), where 
a and 7 are the decay rates of P'{fz) and the experi- 
mental 'P{w) respectively. Since we found in Sec. Ill CI 
that P{F) and P'{Fz) near the bottom are almost iden- 
tical up to a scahng factor {Fz)/{F), the actual decay 
rate of P(/) oc e~"^ is exactly the same as that of the 
(renormalized) P'{fz), so that a = a. Hence, we can 
approximate the exponential decay constant of the force 



distribution as 



To estimate the values of {fic), we worked out two sce- 
narios: we take either p2 — ps or p^ = 0. Together with 
the values of po, pi and 7, taken from the experimen- 
tal data, this yields the values of a listed in the second 
and third column of Tabled Surprisingly, the root mean 
square deviation in a is only 18%, which is rather small 
considering our rather crude estimates of the and the 
fact that Eq. is only approximate. 

Let us briefly recapitulate the discussion above. First, 
we have interpreted the changes in experimental particle- 
wall force distributions of strongly compressed packings 
[m as a change in the packing geometry. To be more 
precise, the overall trends can be understood from the 
expected increase of the number of contacts due to com- 
pression. We demonstrated how one can determine the 
fractions po and pi from the experimental data. Direct 
measurements of these fractions would be very welcome 
as a test of our framework, as well as to extract further 
information of the force distribution P{f). Furthermore, 
our crude estimates in Table H] give reason to believe that 
the force distribution P{f) is actually not much affected 
by the compression. Finally, it seems that for most ex- 
perimental results, where particle deformations are rela- 
tively small, Pq and pi are substantial at the boundary, 
so that Vboundary{w) is similar to Pbuikif) (apart from a 
(5-peak at w = 0). The same argumeiit probably holds 
for recent simulations by Silbert et al. jT^j. 
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2.4 0.23 
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2.29 


0.96 


2.51 
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2.6 0.21 


0.26 


1.60 


1.63 


1.33 


1.96 


37% 


2.8 0.14 


0.18 


1.88 


1.49 


1.54 


1.81 


45% 


3.8 0.00 


0.05 


2.42 


1.57 


1.95 


1.95 



TABLE I: The calculated values for the exponents q, after 
estimating the fractions pn^ from the experimental data of 
Figs. 3a-d of Ref. [lH. The percentage in the first column 
represent the degree of particle deformation. The values of 7 
are taken from Table I of Ref. p4|. 



V. BEYOND THE SIMPLE PICTURE 

In the picture that we have constructed above we char- 
acterize the packing geometry by the fractions , and 
we found that the Vn^iw) are very robust. This is of 



8 



course a vast simplification, since we characterize the lo- 
cal environment of a particle by only one number, namely 
ric- In this section we address the question why this crude 
approach works so remarkably well. For bottom particles 
the situation is particularly simple and insightful, since 
the geometry of the contacts is more or less fixed. There 
is one contact with the bottom, one or two almost hor- 
izontal intralayer contacts and Uc forces from above - 
Fig. 1213. As we have shown in Fig. ^p, the angles of 
these forces display little scatter, so the local texture is 
more or less fixed once ric is given. For bottom particles 
one can thus understand that ric indeed provides a good 
description of the local packing geometry, which justifies 
the decompostion according to Uc. Although for particles 
in the bulk the situation is more complicated, there are 
similar arguments why 7^n^(u') is indeed a robust quan- 
tity, i.e. insensitive for packing geometry. These will be 
discussed in Sec. lV Al We then address the up-down sym- 
metry of the system. Our framework only involves the 
number of contacts from above, ric, and not the number 
of contacts from below, rif,. For bottom particles ric is 
the obvious parameter, but in the bulk of an amorphous 
packing, where the angle distribution is isotropic, there 
is no reason why ric should be more important than nt- 
In Sec. lV BI we therefore investigate weight distributions 
for particles with a given combination {ric, nf,}, which we 
denote by Vn^ntiw). Special attention will be paid to 
particles that have ric 7^ nt, in Sec. IV CI 

A. Why is Vn^ (w) for bulk particles robust? 

It is not a priori clear why Vn^ (w) is rather insensitive 
for the packing geometry, since the definition of Vn^ {w) 
in Eq. (0) involves the joint distribution of the {fi)z that 
push on a particle from above, i.e. P ^(/i)z, • • • , ifnjz^ ■ 
This joint distribution has an explicit geometry depen- 
dence since the projections in the z-direction involve the 
distribution of contact angles ipi . Even if we assume that 
the force magnitude is uncorrelated to its orientation, i.e. 

P[fl,---,fn^) =Pifl,---,fnJH^l,---,^nJ , (10) 

we obtain the distribution of the vertical components 

P {ifi)z, (/njz) by integration over the joint angle 

distribution ^{(fii, ■ ■ ■ , fnj- Therefore, the Vn^iw) have 
an explicit geometry dependence. 

We already saw that this angle distribution is more 
or less fixed for bottom particles. For the polydispersi- 
ties used in this study, the bulk angles have also limited 
room for fluctuations once ric has been specified. For 
example if ric = 3, one typically finds one angle close 
to 7r/2 and two relatively small angles, see Fig. this 
is because the three particles should all touch the up- 
per half of the bead supporting them. Particles with 
ric — 2 also have such an "excluded volume" -like con- 
straint (Fig. 13)), albeit less strong than for ric = 3. Parti- 




90 {|) 180 90 {|) 180 



FIG. 8; (a) For particles with Uc = 3, we plot the probability 
densities for the angles $3(1^1), ^3{ip2) and $3(953), where 
the three angles have been sorted such that ifi < tp2 < V'si 
(b) The probability densities "I>2(vi) and "I>2(</?2) for particles 
with ric = 2; (c) The probability density $i(i/9i) for particles 
with Uc = 1. 



cles with Uc = 1 have an enhanced probability for angles 
around 7r/2, because such contacts make the presence 
of a second contact from above less probable (Fig. 
So, the shape of Vn^ (w) is limited by the geometric con- 
straints on the angle distributions ^{^pi, ■ ■ ■ , <y3n^), which 
are rather peaked. This justifies the picture that the ge- 
ometry dependence of P{w) is mainly due to the pn^ , and 
that the Vn^iw) can be considered invariant. 

Note that the abovementioned constraints on the angle 
distributions imply that the averages {w)na are not sim- 
ply proportional to He- Comparing for example Uc — 1 
and He — 3, we see that the two "extra" forces for Uc = 3 
have a relatively small vertical component; the average 
weight will thus grow less than linearly with Uc- We 
should therefore correct Eq. ^ for the steepness of the 
tails by replacing (ric) with {'uj)n^ ■ Making a cor- 

rection of this type would further refine our analysis of 
the experiment with rubber beads discussed in Sec. lIVBl 



B. Gravity and up-down symmetry 

In our analysis of Viw) we have explicitly broken the 
up-down symmetry, since it only involved the number 
of contacts from above. At the bottom, this is an ob- 
vious choice. Away from the boundary, however, the 
amorphous packings have an isotropic angle distribution 
even though the packings were created under gravity. 
Moreover, we have neglected the term mg in Eq. 
which makes the sum of forces from below equal to the 
sum of forces from above. So in principle one could 
also decompose V{w) according to the number of con- 
tacts from below nt- We therefore investigate Vn^uiiw); 
this can be regarded as a "component" of Vn^iw), since 
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FIG. 9: (a) Vi2{w) (solid line) and P2i(w) (dotted line); 
(b) Visiw), P22{w) and Vniw); (c) T23{w) and V32{w); (d) 

Fig-Eb shows that Visiw), V22{w) and ■P3i(w) are al- 
most identical. The same holds for V2z{w) and 'Ps2{w) 
(Fig. 15];) , so the total coordination number ric + rn, ap- 
pears to be a more fundamental quantity than just Uc or 
rif,. Fig. 131 furthermore shows that the quadratic scaling 
of 'P33(w) is somewhat more pronounced than for 7^23 (w^) 
and V-i2{'w)\ it seems that the presence of 2 contacts from 
above or below inhibits the pure quadratic scaling. 

The presence of gravity is noticed, however, for 7^12(10) 
and V2i{w) which do show some differences (Fig. IHt)- 
These particles have only 3 contacts and were less re- 
stricted during the formation of the static force network 
by the "cage" surrounding them. This allowed gravity to 
influence their final movements more than for particles 
with nc + Tih > 3. Obviously, this effect is even stronger 
for particles with only 2 contacts, which typically have 
K,nb} = {0,2}. 

To further investigate the up-down symmetry, we list 
the fractions Pn^n,, of particles with a certain ric and n;, in 
Table Hll For all particles with 3 or more contacts these 
fractions are almost perfectly symmetric. From this we 
conclude that in the amorphous packings, the up-down 
asymmetry due to gravity is only noticed by particles 
that have 2 or 3 contacts. 



C. Particles with Tie 7^ ni, 

We have seen that for particles with {ric, nf,} = {3, 1} 
or vice versa, the small weight behavior is ~ w^, which 
is different from the scaling predicted by Eq. ©. This 
breakdown of our simple picture can be understood as 
follows. A particle that has 4 contacts can either have 
{uc^nb} = {3,1}, {nc,nb} = {2,2} or {ric.nb} = {1,3} 
depending on the precise orientations of the forces with 
respect to gravity. However, if we were to define the 
weights by projecting the Fij at a small angle with re- 
spect to gravity, a particle with 4 contacts can easily 
change from {ric, n^} = {3, 1} to {2, 2} or even to {1, 3}. 
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TABLE II: Fractions pn^n^ expressed in percentages; the 
numbers are almost up-down symmetric, except for rattlers 
(particles with 2 contacts). From these fractions one finds the 
average coordination number {ric + rib) = 4.51. 



However, we have seen that there is no "preferred" pro- 
jection direction, since gravity has only very little effect 
on our packings. Hence, it is not surprising that the 
'Pubiia (^) depend on 4- n\, and not on or rib individ- 
ually. 

But what determines the precise scaling for small 
weights? Consider a particle i with ric = 3 and = 1. 
The three forces pushing it from above, Fn^ Fi2 and F^^ 
are not independent: force equilibrium in the direction 
perpendicular to Fa (the force pushing from below) re- 
quires {Fh + Fi2 + Fi3^ • n = 0, where Fi4 ■ n ^ 0. This 

reduces the number of independent forces from above to 
only 2, since the third is determined by mechanical equi- 
librium. As a consequence, the scaling behavior for small 
w will be 7^31 (w) oc w. 

For particles with = 3 and rib — 2, the 5 forces 
are also coupled through mechanical equilibrium. In 
this case, however, one can not distil a relation between 
the forces from above only, such as we did for parti- 
cles with {ric, rib} = {3,1}. So one still expects that 
7^32 (w) oc w^, as is observed in Fig. |^. Nevertheless, 
this illustrates that two dimensional mechanical equilib- 
rium does introduce correlations between all forces push- 
ing from above. This limits the validity of our argu- 
ments used in Sec. IIIII for bulk particles. At the bottom 
our analysis is still valid: horizontal equilibrium can be 
accomplished by the forces between neighboring bottom 
particles (see Fig.Eb), so the forces from above can really 
be considered as independent. 

D. Summary 

In this section we have addressed the limitations of 
our simple geometrical framework. We have shown that 
the observation that Vn^ {w) is insensitive to packing ge- 
ometry originates from excluded volume-like correlations 
between the angles at which forces press upon a bead 
(Fig. IS)) . This is the subtle underlying reason why our 
simple picture, where we characterize the local packing 
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geometry by only one number Uc, is good enough to in- 
terpret experimental and numerical data. We have fur- 
thermore studied the effect of gravity by decomposing the 
weight distribution according to the number of particles 
from below {rib) as well. We found that gravity breaks 
the up-down symmetry only mildly in our simulations; 
the distributions Vn^n^ (w) depend on the coordination 
number ric + Ub rather than on Uc or Ub independently 
(Fig. I^l . This further refines the analysis of the relation 
between packing geometry and force network statistics in 
the bulk of a packing; at the boundary, it is sufficient to 
consider only the number of contacts from above (ric). 



of Fij. Since these interparticle forces are more funda- 
mental than the weights, we investigate the statistics of 
the quantity qW in Sec. IVI Al this will shed new light 
on the discrepancy for small forces between the q-model 
and experimental data. In the light of our finding that 
the contact geometry and in particular ric plays a cru- 
cial role, the standard g-model is clearly limited since it 
fixes Uc- In Sec. IVI B1 we therefore extend the g- model to 
have randomness in its connectivity (i.e. to allow for a 
range of ric's), and find that, as expected, the V{w) can 
be manipulated by changes in the connectivity. 



VI. WEIGHT AND FORCE DISTRIBUTIONS IN 
THE g-MODEL: THE ROLE OF CONNECTIVITY. 

In this section, we investigate to what extent the re- 
sults obtained for the Hertzian sphere packings can be 
understood within the context of the (/-model and its 
generalizations. In the standard version of the model, the 
particles are positioned on a regular lattice, and the par- 
ticle weights are stochastically transmitted to the neigh- 
bors in the layer below [5[. The weight on a particle i 
splits up into Uc fractions qij, and the total weight ex- 
erted on a particle j in the layer below then becomes 

ly, =mff + ^(7yW„ (11) 

i 

where the term mg can be neglected at large depth. The 
fractions qij obey the constraint 

E9u = l' (12) 
j 

which assures mechanical equilibrium in the vertical di- 
rection. They can in principle also be deduced for 
more realistic packings: from definition one finds 

The simple form of the g-model has allowed for a num- 
ber of exact results of which the most important is the 
solution for the uniform (/-distribution. This uniform q- 
distribution assigns an equal probability to each set of 
{qij} that obeys Eq. (|12|l . and serves as a generic case. 
The rescaled weights w then become distributed as 

7'„,(H = cu;"=-ie-"=™, (13) 

where ric is fixed for a given lattice, and c is a normaliza- 
tion constant. Note that these solutions have the same 
qualitative behavior as those found in our molecular dy- 
namics simulations: for small weights Vn^iw) oc w"''"^, 
and the probability for large weights decays exponen- 
tially. 

The (/-model is thus an effective minimal model for 
the weights W. It is clear that the product of qij and 
Wi has a natural interpretion as the vertical component 



A. Distribution of interparticle forces: P{q'w) 

A direct comparison of Eqs. ^ and (|ll|l shows that 
the product qijWi has a natural interpretation as the ver- 
tical component of fij . Since the interparticle forces are 
more important than the weights, it is interesting to in- 
vestigate the statistical properties of the bond quantity 
qw. To obtain the distribution P{qw), let us start with 
the transformation from P(qw) to Vn^iw): 

VnA^) = ^ d[qw)iP l^{qw)i) ■ ■ ■ 

X rf(gw)„^p((gu))„J 

xs(w~Y.{qw)}j. (14) 

Here we assumed that the {q'w)i are uncorrelated, which 
is valid for the uniform g-distribution [l^. For the cor- 
responding Laplace transforms, denoted by P{s) and 
"Pn^is) respectively, this relation becomes 

-PnAs) ^ {P{s)y^ ■ (15) 

Since the Laplace transform of Eq. H13|l is of the form 
1/(1-1- 3)""=, the distribution of qw reads: 

Pis) = => P{qw) = e-«'". (16) 

We thus find (for the uniform (/-distribution) that P{qw) 
is a pure exponential, independent of the number of con- 
tacts lie- Again, this is very similar to the results for our 
Hertzian sphere packings: the distribution of "interpar- 
ticle forces" P{qw) is finite for small forces, whereas the 
distribution of weights depends on ric as given by Eq. 
Moreover, this resolves the discrepancy for small forces 
mentioned in the introduction: the q-model predicts a 
vanishing probability densitity for small weights, but not 
for small forces. 
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FIG. 10: The g-model with a random connectivity: (a) with 
a probabihty p we cut one of the three bonds; (b) bottom 
effect in the g-model with random connectivity. In the bulk 
p = 0.3 (open circles) and at the bottom p = 0.9 (dots); this 
corresponds to {po, Pi, P2, ps} ~ {0.00,0.03,0.24,0.73} and 
{0.03, 0.19, 0.44, 0.34} respectively. 



B. Including geometry effects 

From Sec. IIIII it is clear that the weight distribution 
V{w) in Hertzian sphere packings is very sensitive to the 
local packing geometry. Since the g-model is defined on a 
regular lattice, with fixed connectivity, it can not capture 
the behavior of V{w) in disordered packings with fluctu- 
ating ric- This extra degree of disorder can be included, 
for example, by "cutting" some of the bonds of the reg- 
ular lattice. We illustrate this with the 2-dimensional 
square lattice depicted in Fig. llOb . For each site, the 
weight is transmitted downwards through either 2 or 3 
bonds with probabilities p and 1 — p respectively; in the 
former case we randomly cut one of the available bonds 
and generate the two remaining qij according to a uni- 
form distribution satisfying Eq. (|12|l . This generates par- 
ticles with Tic — 0, 1, 2 and 3, since all bonds arriving at 
a site have a probability of p/3 to be missing. For sim- 
plicity, we introduced the disorder in ric by means of one 
parameter p only; as a consequence, we can only obtain 
a limited set of {pn^}- 

With this model, we have tried to mimic the bulk- 
bottom behavior of 7^(?x;) that was observed in the amor- 
phous packings (Fig.|3Ji). In the bulk layers we took out 
bonds with probability p ~ 0.3, and for the bottom layer 
we took p — 0.9; the result is shown in Fig. llOb. Indeed, 
the change in the fractions is sufficient to reproduce a 
transition of V{w) reminiscent of what has been observed 
in our Hertzian sphere packings (compare with Fig. [S^). 

C. Conclusions for the g-model 

Although it is known that the g-model does not prop- 
erly describe the spatial structure of the force-network 
[l2|. it remains a very instructive theoretical framework 
for the statistics of force fluctuations. While in the stan- 
dard case the disorder in the system is represented by the 
stochastic fractions only, we have shown that when 
also the connectedness is chosen to be random, the model 
displays most features of realistic packings. 



Let us conclude this section by mentioning that the 
idea to leave out some of the bonds of a regular lattice 
is not new 20]. In these studies, however, bonds were 
cut in a particular manner to build up directed force 
chains. We have shown that such long-ranged structures 
are not important for the behavior of V{w), since they 
only depend on the local packing geometry. 



VII. TOP-DOWN RELAXATION OF 
FLUCTUATIONS 

So far, the discussion has been limited to situations 
well below the top surface of the packings. The data 
of the Hertzian sphere simulations were taken at least 15 
layers below the top surface and the results of the g-model 
(presented in the previous section) all correspond to the 
limit of large depths. In this section we investigate the 
top-down relaxation of the force and weight distributions. 
At the top surface of the Hertzian sphere packings, there 
are only weight fluctuations due to grain polydispersity. 
The question we address is how fast the force and weight 
fluctuations build up towards a bulk distribution, as a 
function of depth. 

These results can then be compared to the relaxation 
in the g-model. Interpreting the downward direction as 
time, this corresponds to transient behavior towards the 
"stationary" solutions given in Eqs. H13|l and H16|) . This 
top-down relaxation of fluctuations forms an additional 
test to qualify various theoretical models, very much like 
the Green's function measuring the response to a local- 
ized load on the top-surface In our case, we start 
from spatially (nearly) homogeneous conditions in the 
top layer and see how fluctuations build up. 



A. Top-down relaxation in Hertzian sphere 
packings 

A good way to quantify changes in 'P{w) and P{f) is 
to study their second moments (w^) and (/^). For a dis- 
tribution of zero width these second moments are unity, 
and they increase as the fluctuations become larger. In 
Fig. ^2 we show the second moments as a function of 
the height h, which is deflned as the distance from the 
bottom boundary. Since the packings are strongly dis- 
ordered, the precise location of the top surface will be 
slightly different for each realization; it turns out to be 
located around h = A6. 

Let us flrst consider the broadening of the weight dis- 
tribution shown in Fig. lllb . As already mentioned above, 
the weight fluctuations at the top surface are entirely due 
to polydispersity of the grains. Using a flat distribution 
between 0.4 < r < 0.6 this corresponds to (w^) « 1.11, 
which is consistent with our simulation data. The sec- 
ond moment approaches its bulk value already at a depth 
of approximately 10 particle diameters. The figure also 
shows the sharp transition of V{w) at the bottom bound- 
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FIG. 11: The second moments (a) {w^) and (b) (/^) as a func- 
tion of height h in simulations of Hertzian sphere packings. 
The arrow indicates the location of the top surface, around 
h = 46. Both for the forces and the weights one finds a fast 
top-down relaxation of the moments. 



ary. The second moments of P{f) are shown in Fig.lllb. 
One again observes a relaxation over approximately 10 
layers, towards a bulk value; P{f) does not change sig- 
nificantly near the bottom boundary. Note that both the 
force and weight distributions become slightly narrower 
as the depth increases below heigths of the order of 30. 
This may be attributed to an increase in particle defor- 
mations [l5l| . 

We thus find that the typical length scale for force and 
weight fluctuations to saturate is approximately 10 parti- 
cle diameters. This provides another important criterion 
to distinguish between difi'erent theoretical models. 



B. Top-down relaxation in the q-model 

The top-down relaxation is well understood for the q- 
model without the so-called injection term, i.e. mg = 

m Eq. ini) Ellis. Before extending these results to 

the g-model with injection, we briefly recapitulate the 
results of the q-model without the injection term mg. 
This version of the model can be interpreted as a pack- 
ing of weightless particles, supporting a homogeneously 
applied force. To distinguish between the g-model with- 
out injection from the model with injection, we denote 
the weight distributions at depth t by TZ^^^w) (without 
injection) and by (with injection). 

For the uniform (/-distribution, it has been shown that 

M 



d-l 



for t 



oo, 



(17) 

where d is the dimensionality of the packing. The sta- 
tionary solution V(w) is given by Eq. p3|l and F{w) is 
the shape of a typical deviation. It is clear that all second 
and higher order moments {w'^) approach their asymp- 
totic values according to the same power-law. This slow 
relaxation towards V{w) is caused by diffusion of cor- 
relations, which takes place in the {d — l)-dimensional 
correlation space p^ . 



Let us now investigate how the injection term mg af- 
fects the top-down relaxation. We first note that the 
recursive relation for the weights, Eq. I|ll(l . is a linear 
equation. The q-model with injection can therefore be 
interpreted as a superposition of g-models without injec- 
tion, with differently positioned initial layers. Although 
it is not a priori clear how this superposition property 
is refiected in the weight distributions V^^^w) (with in- 
jection) and TZ^*\w) (without injection), we propose the 
following approximate mapping: 



^ ^ t+1 ^ ^ 



(18) 



If we combine this with the exact result of Eq. (|17|l . we 
obtain the following relaxation as i — > oo: 



-= d = 2 



log(t) 



t 



d = 3 



d> A 



(19) 



This relaxation behavior is indeed observed in our numer- 
ical simulations with d = 2 and d — 3, using a uniform 
q-distribution. In Fig. ^1 we show the results for an fee 
packing {d = 3). We plot t x \{w'^)^*^ - 4/3| as function 
of depth t, where (z/;^)(°°^ — 4/3. The climbing straight 
line on the lin-log plot confirms the remarkable \og(t)/t 
relaxation. We also plot the same data for the g-model 
without injection; this curve becomes fiat in agreement 
with Eq. ((T7jl . 

Although the mapping of Eq. (|18|l is definitely not ex- 
act, it apparently captures the main physics of the relax- 
ation process. This can be understood as follows. There 
are two slow processes involved: (i) the increasing num- 
ber of layers reduces the contribution of each layer of 
"injected" weights effectively as l/t; (ii) each layer of in- 
jected weights relaxes as {1 / individually. Nat- 
urally, the total relaxation is dominated by the slower of 
these two processes. In the special case of d = 3 both 
powers are l/t, leading to a logarithmic correction. Fi- 
nally note that since the downward g-values are statis- 
tically independent from the weights, the "force" fiuctu- 
ations simply follow from {[qwY) = {q^){vP'), and thus 
display the same relaxation as the weights fiuctuations. 



C. Conclusions concerning top-down relaxation 

We have studied the top-down relaxation of the sec- 
ond moments {uP') and (/^), which quantifies how 'fast' 
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t|<w2> -4/3 I 




1 10 100 1000 

t 

FIG. 12: Relaxation of the second moments with injection 
(climbing line) and without injection (flat line) towards their 
asymptotic values 4/3 in the 3D 5- model. Since we plot t x 
I (w^)^*'— 4/3| along the vertical axis, the climbing straight line 
confirms the log(f)/t relaxation for the g-model with injection. 
Without injection the relaxation is simply 1/t. 



the weight and force distributions approach their bulk 
shapes. The g-model predicts a power-law relaxation 
with a logarithmic correction for 3D packings, Eq. (|19|l . 
However, we find no evidence for such a sfow relaxation 
in our simulations of Hertzian spheres, which indicate 
that a bulk distribution is reached after approximately 
10 layers of particles (Fig. Ill|l. In the g-model with in- 
jection, for example, the second moment after 10 layers 
still differs around 20% from its asymptotic value. 

Let us provide two possible explanations why the q- 
model fails to describe this relaxation process. A first 
problem of the model is that it assumes some fixed q- 
distribution 77(g): we have seen that the q's can in prin- 
ciple be derived from the forces as = (^y) /^ii so 

a relaxation in P{f) and Vlw) should result into a re- 
laxation of ri{q) itself. This clearly shows the difficulty of 
encoding the force behavior into a stochastic variable q in 
a self-consistent manner. Another problem of the model 
is that it assumes a top-down propagation of forces. The 
up-down symmetry is therefore broken explicitly in the 
g-model, whereas in our Hertzian sphere packings we find 
only a very weak symmetry breaking. In principle, force 
networks are defined by the equations of mechanical equi- 
librium, which generically are underdetermined |23l |24| 
and hence can not be solved by an iterative (top-down) 
procedure. Instead, one has to solve this set of coupled 
equations "simultaneously" for all particles in the sys- 
tem, and except for the (small) mg term, there is a nat- 
ural up-down symmetry in this system. The absence of 
this up-down symmetry in the g-model could of course 
strongly affect the top-down relaxation. 



VIII. DISCUSSION 

We have shown that in order to understand the statis- 
tics of force networks, it is crucial to distinguish between 
forces and weights. We have found in our simulations 
that the force distribution P{f) is very robust, in the 
sense that its shape does not depend on details of pack- 
ing geometry. The weight distribution V{w), on the other 
hand, is very sensitive for the local packing geometry. We 
have demonstrated that a decomposition according to the 
number of contacts that press on a particle from above. 
Tic, is sufficient to understand this geometry dependence. 
Reinterpretin g e xperiments on strongly deformed rub- 
ber particles [ll| within this framework, we find strong 
evidence that P{f) essentially remains unaffected even 
by very large particle deformations. To further test our 
framework experimentally, one can manipulate the num- 
ber of contacts at the boundary by placing a layer of 
relatively small or large beads at the bottom. For small 
beads, the fractions po and pi will be enhanced, leading 
to a large V(w) for small w, and a slow exponential decay 
for large w. Relatively large bottom beads should lead 
to a V{w) that is strongly peaked. 

The present work provokes a number of questions. 
First, we observe that most of our simulation results, 
like the shapes of P'{fz) and V{w)^ can be understood in 
terms of local packing geometry only. This suggests that 
long-range correlations are not dominant, at least not for 
the 'one point' force, weight and angle probability distri- 
butions. We therefore question whether the behavior of 
P{f) observed at the jamming transition |^ Q reflects 
a long-range structural change of the force network. In 
particular, we expect that the role of 'force chains' can 
only be understood from two or more-point correlation 
functions, and not from P{f) only. 

A related problem is that the q-model fails to describe 
problems that involve spatial structure of the force net- 
work. Although the model is able to capture many fea- 
tures of force and weight statistics (Sec. IVI|I . it does not 
produce the top-down relaxation oiV{w) that is observed 
in the more realistic Hertzian packings. Alongside with 
the incorrect prediction of the response function [T^ , this 
indicates that spatial dependence is not correctly incor- 
porated within the g-model. This may be due to the fact 
that, in general, recursive models do not acknowledge the 
structure of the equations describing mechanical equilib- 
rium. These equations are typically underdetermined 2"^ 
and cannot be solved in a recursive manner. In a recent 
paper 0|, we therefore propose a different theoretical 
approach, in which we start from the equations of me- 
chanical stability and exploit the undetermined degrees 
of freedom. 

Another important issue for future study is clearly the 
role of friction and dimensionality. Our numerical study 
has been done in two dimensions with frictionless spheres; 
however, recent studies indicate that the coordina- 
tion number for 3D packings with friction is similar to 
those of 2D frictionless packings. Qualitatively, the pic- 
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ture we have advanced is therefore expected to capture 
the reahstic case of three diraensions with friction, be- 
cause our phase space arguments are independent of di- 
mension. 
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APPENDIX A: LOGARITHMIC DIVERGENCE 

OF P'if.) 

In Sec. Ill CI we encounter the following integral: 
1 



dip- I dfPif)dif,-fsm{^)) 
Jo 

2 1 



The function P{f) represents the probability density 
function of / = |/|, which we can assume to be regular 
on the entire interval (see Fig. |2) . The behavior for small 
fz is not trivial, since the integrand diverges at the lower 
bound of the integration interval. For each non-zero fz 
this does not lead to a singularity, since 



P'Uz) 



2 df 
2 



P{f) 



TT 



du- 



Pjufz) 



(A2) 



The integral over l/\/v? — 1 is convergent for u ^ \ and 
the function P{ufz) falls of fast enough as (w/z) oo. 
For fz = 0, however, the integral diverges as u — > oo. To 
obtain the asymptotic behavior we rewrite the integral 
as 



2 f° 

P'U.) = - / 



du 



P{ufz) 



2 f°° /I 1 

- / duPiufz) -== - 
Jl V V M - 1 " 



.(A3) 



The second term is convergent since the term between 
brackets behaves as in the limit m ^ oo. We thus 
find that 

p'ifz^o) ^ - r df^+o{\) 
^ jf. f 

--P(0)ln(/z) + O(l) . (A4) 

TT 



APPENDIX B: RELATION BETWEEN TAILS OF 

P'(/,) AND PnAw) 

In this appendix we derive the large weight behavior of 
Vn^ (w) from the tail of P'{Fz), assuming that the various 

(^Fij in Eq. © are uncorrelated. We consider decays 

both faster and slower than exponential, of the form 



P'{Fz) oc e-"-^-^/^^^^'' for i^^ ^ oo . (Bl) 
We show that, after rescaling (w) to unity, this leads to 

VnAw)oce-^^\ (B2) 

with 



7 




(B3) 



This means that the tail of the weight distribution is of 
the same nature as that of the forces, but with a different 
prefactor 7. The tails get steeper for increasing ric, since 
the reduced probability for small w (due to a lack of phase 
space) must be compensated to keep (w) = 1. 

The above results are obtained as follows. Rescaling 
all forces in Eq. ^ a,s Xi — {Fz)^/W, one obtains the 
probability for large weights 



dxi ■ ■ ■ dxn^ e <-p^> 



(B4) 

where S denotes the hyperplane 1 — with all Xi>Q. 

For /? > 1, the probability density on S has a maximum 
at Xi = ^jric-, which becomes sharply peaked for increas- 
ing W . Physically, this means that the dominant contri- 
bution for large weights will come from all Fz being equal, 
namely W juc- Approximating the integrand by a Gaus- 
sian around its maximum value, we find that the "width" 
decreases as a power of W only, namely l/iy("<:~i)'^/2. 
Hence the leading behavior for large W is given by the 
maximum value of the integrand, i.e. e T^^'^ /("<=) 

For /3 < 1, the probability density has a minimum at 
Xi — Xjuci and the dominant contribution now comes 
from Xi = \ and Xj^i = 0. This means that typically 
only one of the forces accounts for the whole weight. The 
part of the integral around Xi = \ can be approximated 
by 



dxi ■ ■ ■ dxji e 



(B5) 



where Se denotes the part of S for which 1 — Xi < e. 
This approximation becomes exact for — s- 00 as long 
as W^e < 1; we take e = l/W^'^ with < S < 1 - (3. 
Working out the integration over S^, one finds 
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e 



dye <^-> 



(B6) 



as oo. The part of the integral outside the areas 

Sf is smaher than W"''~^e TpIW'^ (i+W' ) t\ms 



be neglected. So also for /3 < 1, the leading behavior 
for large W is simply given by the maximum value, i.e. 

As mentioned in Sec. IIIII the Vn^iW) obtained by 
Eq. |SJ) are not properly normalized, since (W) = {fz)nc. 
If we rescale the average weight to unity, we obtain the 
results of Eqs. and (|B3)| . 
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